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ABSTRACT 

The stalling radius of a merging massive binary black hole (BBH) is expected to be below 
0''1 even in nearby galaxies (Yu 2002), and thus BBHs are not expected to be spatially re- 
solved in the near future. However, as we show below, a BBH may be detectable through the 
significantly anisotropic stellar velocity distribution it produces on scales 5-10 times larger 
than the binary separation. We calculate the velocity distribution of stable orbits near a BBH 
by solving the restricted three body problem for a BBH embedded in a bulge potential. We 
present high resolution maps of the projected velocity distribution moments, based on snap- 
shots of ^ 10^ stable orbits. The kinematic signature of a BBH in the average velocity maps 
is a counter rotating torus of stars outside the BBH Hill spheres. The velocity dispersion maps 
reveal a dip in the inner region, and an excess of 20-40% further out, compared to a single BH 
of the same total mass. More pronounced signatures are seen in the third and fourth Gauss- 
Hermite velocity moments maps. The detection of these signatures may indicate the presence 
of a BBH currently, or at some earlier time, which depends on the rate of velocity phase space 
mixing following the BBH merger 

Key words: black hole physics - galaxies ; nuclei - stellar dynamics. 



1 INTRODUCTION 

The discovery that most galaxies harbour a massive black hole 
(BH) at their core (Magorrian et al. 1998), and the commonly ac- 
cepted interpretation of cosmological structure formation simula- 
tions, that galaxies grow by mergers (e.g. Kauffmann et al. 1993, 
cf. Dekel & Bimboim 2006 that growth is mostly by gas accre- 
tion) implies that massive binary black holes should commonly 
form at the core of galaxies (Volonteri et al. 2003). Unlike the pre- 
vious view, that the formation of such a binary results from the 
rare process of merging during the life of the rare objects, quasars 
(Begelman, Blandford & Rees 1980). Dynamical friction on the 
background stars leads to inspiral of the binary on a dynamical 
time-scale, until the binary becomes hard, which occurs at a binary 
separation close to 
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where G is the gravitational constant, M. is the mass of the pri- 
mary BH, q is the mass ratio of the BHs (g ^ 1), and a is the bulge 
ID velocity dispersion. The inspiral rate slows drastically when 
'' ^ '^h (6-g- Merritt 2006a), and further inspiral is set by the rate at 
which stars diffuse in phase space into the loss cone (Frank & Rees 
1976)' . The minimum phase space diffusion rate is set by the two 
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as pointed out by 



body relaxation rate. This mechanism may be fast enough in the 
lowest luminosity bulges, where the stellar cores are the densest, 
and may lead to merger through gravitational wave emission in less 
than the Hubble time (but see recent suggestions that the core den- 
sity in the Milky Way is smaller than originally thought, Merritt 
2010 and references therein) . In more luminous bulges, the binary 
inspiral is expected to stall (e.g. Makino & Funato 2004; Merritt 
2006b, and references therein), leading to the final parsec prob- 
lem (Begelman et al. 1980). The largest plausible stalling radius in 
the nearest galaxies is just below Of 1 (Yu 2002), and is thus gen- 
erally unresolved. The stalling may be overcome by diffusion of 
orbits into the loss cone induced by tangential forces. Either due to 
bulge triaxiality or bar like structure, or by massive perturbers in 
the form of molecular clouds or a third BH due to another merger 
(e.g. Merritt & Poon 2004; Berczik et al. 2006; Hoffman & Loeb 
2007; Alexander 2007; Perets & Alexander 2008). Alternatively, 
inspiral may be induced by angular momentum loss of the BBH 
to circumbinary gas (e.g. Armitage & Natarajan 2005; Mayer et al. 
2007; Cuadra et al. 2009). These processes may lead to a merger 
on time-scales well below the Hubble time. 

The formation process of the BBH ejects stars from the bulge. 
This occurs on scales significantly larger than a^^, and it flattens 
the stellar density profile. This BBH "scouring" may be respon- 
sible for the formation of the core present in luminous ellipticals 
(Milosavljevic & Merritt 2001; Merritt & Milosavljevic 2005). Re- 
cent high quality observations by Kormendy & Bender (2009) indi- 
cate a remarkably tight congelation between the deduced mass defi- 
ciency in the core region, and the black hole mass. This correlation 
is likely a signature of the cumulative scouring effects of merging 
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BHs during the BH growth. Are there any other signatures induced 
by the merger? 

If the binary stalls, then a significant fraction of galaxies may 
harbour a compact BBH (e.g. Volonteri et al. 2003), which cannot 
currently be directly resolved (Yu 2002). The binary affects the 
stability of stellar orbits on scales up to several times larger than 
the binary separation. The presence of a binary may thus be in- 
ferred through its signature on the background stellar distribution 
and kinematics (e.g. Cappellari & McDermid 2005; Kandrup et al. 
2003). The purpose of this paper is to determine this signature and 
provide a tool which can be used to detect the effect of a binary 
black hole on scales 5-10 times larger than aj,. ^ scale which may 
be resolved in nearby galaxies with current angular resolutions. 

We first explore the stability of orbits near a BBH by calcu- 
lating stability maps. These maps present the time it takes for a 
test particle to become unbound, as a function of its position in 
velocity space, at a given launching radius. A similar concept ap- 
pears in Wiegert & Holman (1997) in the context of small bodies 
near a binary stellar system, but the coordinates used there were the 
semi-major axis and the inclination of the orbit, which are less rel- 
evant for the kinematic signature explored here. The stability maps 
elucidate how the velocity phase space populated by stable orbits 
evolves with integration time, and how it varies with distance from 
the binary. The maps shows how the velocity distribution function 
/(v) evolves from a smooth isotropic function at a distance r^ a\^, 
to a highly anisotropic function, which produces the BBH kine- 
matic signature. 

We derive the BBH kinematic signature by integrating orbits 
of test particles around a massive BBH, i.e. by solving the restricted 
3-body problem, for a binary in a circular orbit, embedded in a 
bulge potential. We take snapshots of the system once it approaches 
a steady state, and produce maps of the projected distribution of the 
stars and their velocity distributions. Similar scattering experiments 
were already done in the past (e.g. Hills 1983; Mikkola & Valtonen 
1992; Quinlan 1996; Sesana et al. 2007). The purpose of these stud- 
ies was to derive the effect of the stellar background on the BBH 
merger rate, while our purpose is the reverse, to derive the effects 
of the BBH on the background stellar population. The latter calcu- 
lations were already carried out by Milosavljevic & Merritt (2001) 
based on A'-body simulations. These simulations were limited to 
N <\Q and cover a region more than 10 larger than a\^, thus 
they do not allow accurate mapping of the projected line of sight 
velocity distributions (LOSVDs) on scales of a few times flj,, as 
done here. Since the background stellar mass within flh is neg- 
ligible (see below), and 2-body scattering cross sections are still 
small, A'-body simulations can be replaced by the much faster scat- 
tering experiments. These experiments allow us to probe the stellar 
kinematics with a roughly 10 higher resolution, in the relevant re- 
gion, compared to Milosavljevic & Merritt (2001). In §3 we present 
the method of the calculation, and in §4 we present results on the 
phase space stability regions, maps of the projected velocity distri- 
bution moments, velocity moments along some slit positions, vari- 
ous LOSVDs, and also stellar density profiles and projected density 
maps. The results are discussed in §5, and summarized in §6. 



2 MODEL 

2.1 Orbital Elements 

The two BHs are assumed to be in circular orbits in the x-y plane, 
about their centre of mass at (x,y) = (0,0). The orbital period, the 



relative velocities, and radii are: 
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where a is the binary separation, corresponding to the stalling ra- 
dius. We take a = a^^, i.e. the stalling separation is the hard bi- 
nary separation (defined in equation 1). The recent simulations 
of Merritt (2006a) (table 1 there) indicate that the stalling radius 
agrees with a^ to typically 20%. Combining equations (1) and (3) 
gives: 



V _ 2{i+q) 
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The physical scales are set by specifying two parameters, 
say M. and V. However, since M, and (T are correlated through 
the M-a relation (Gebhardt et al. 2000; Ferrarese & Merritt 2000; 
Tremaine et al. 2002, and citations thereafter), only one parameter 
is required. We use the recent Giiltekin et al. (2009) relation: 

logM8=0.12±0.08 + (4.24±0.41)log(T200 (7) 

where Mg =M./10'* M;, and a20() = (7/200 km s"', which gives: 

ah = (3.1±0.3)^M"«±"05pc (8) 

= (3.5±0.3)^cT|fo4±o.'»pc. (9) 

Thus, the binary separation is expected to be on the parsec scale. 
Similarly, the orbital period is: 



T = (5.5±0.3)xl04-^^M«-2''±0'«yr 



(1+?) 



(5-0±0.3)xlO'*^a]or^»'yr. 



(10) 
(11) 



Thus, a characteristic orbital time-scale of lO'* years, with a rela- 
tively weak dependence on mass. 



2.2 Model Units 

For the sake of simplicity, we use the following values G = 1, 
M, = 1, and a = 2. To convert back to physical units, one needs 
to specify the physical value of M. ; and then the length, which is 
measured in units of a^/Z, is given in parsecs by equation (8). Time 
is then measured in units of (aj^/8GM.)' " and velocity in units of 
(2GM,/ah)'"- Alternatively, the conversion from dimensionless 
to physical units can also be made using the value of (T2oo- 

Note that models with different q correspond to different phys- 
ical scales for the same mass scaling (see equation 1). Table 1 gives 
the physical scales of the model units for M. = 10 M©, and q=\ 
and 0. 1 used below. Also note that a in model units is a function of 
q only (equation 6). 
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Table 1. The orbital elements in model units and physical units. The con- 
version is for M, = 10^ M,:,, and the conversion factors are listed in the 
lower part of the table. 





Expression 


Model 


= 1 
Physical 


Model 


= 0.1 
Physical 


a 
T 

V 




2 
12.57 

1 

0.25 


1.55 pc 

12 783yr 

745 km/s 
186 km/s 


2 
16.94 

0.74 
0.11 


0.28 pc 
1 336 yr 

1 296 km/s 




/ '/ 


186 km/s 




V 8(1+'/) 




Length 

Time 

Velocity 


1 

1 
1 


0.78 pc 

1017 yr 

745 km/s 


1 

1 
1 


0.14 pc 

79 yr 

1 747 km/s 



2.3 Bulge Properties 

The BBH is embedded in an isothermal sphere, i.e. the veloc- 
ity distribution of the stars is a Maxwell-Boltzmann distribution 
/(v) = Avexp^'' /-'^", where v is the magnitude of the velocity, 
A is a normalization coefficient. The assumed stellar density pro- 
file is p{r) = cy- /iKGr-, the self-consistent solution for a singular 
isothermal sphere. This is only an approximate solution for a finite 
mass bulge, and is not appropriate within the BH sphere of influ- 
ence. Here we use /(v) and p{r) as convenient initial conditions 
for the stellar distribution close to the BBH. Also, p{r) allows a 
convenient representation of the bulge potential, and is surprisingly 
accurate for the average properties of massive (> 3 x 10 Lp,) el- 
lipticals, within their effective radius (Koopmans et al. 2009). To 
avoid the non physical divergence of p(r) at r = 0, we assume a 
core structure: 
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where h is an arbitrary break radius and 

The integrated bulge mass derived from the above density field is: 
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The expression for the gravitational potential (i.e. the bulge poten- 
tial) is: 
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The /V-body merger simulations of Milosavljevic & Merritt 
(2001) indicate that at the time the binary becomes hard, for a short 
while, an r^^ density profile extends down to the scale of ajj^- We 
therefore assume h = a/2 (in model units, h= I, but note the final 
larger core radius produced at the end of the simulation). From here 
on, we work only in model units (defined in §2.2). 

The calculations of the binary orbital elements (equations 2-5) 



Though this is not necessarily a realistic result, given the small number 
of 10 stars interacting with the BH in their simulation within ri„fl, and the 



ignores the bulge mass within the BBH orbits. The enclosed mass 
within the radius of the secondary BH (equation 5), which resides 
outside the uniform density core, is: 
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Thus, M(i?2) ^ 0.04M,, and the bulge mass was therefore ne- 
glected in the derivation of V and T made above. 

The bulge mass grows linearly with r, and becomes larger than 
the combined BH mass at: 
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implied very short relaxation time. 



This definition of the radius of influence differs somewhat from 
the commonly used definition of G/W./cT", or 8 /g + 8 in the model 
units. For ^ = 1, rjng ~ 17 and for q = Q.\, rjnfl ~ 49; but note that 
in physical units the latter number is smaller (see Table 1). The 
apparent divergence of rjnfl for q^ Q results from the divergence 
of K as Oh ^ (equations 1, 3). Below we also simulate orbits 
around a single BH, for the purpose of comparison. This simulation 
is designated as the q = Q case, but we do not use the V based 
normalization, due to its divergence. 



3 METHODS 

3.1 Orbit Integration 

The orbit of each test particle (representing a star) is solved using a 
5 order Runge-Kutta method with adaptive step size control. The 
fractional error tolerance in the code was set to 10^ , further lower- 
ing this value had no detectable effect on the results. The accuracy 
of the calculation was also verified through conservation tests of 
the value of the Jacobi integral, a constant of motion in the circular 
restricted 3-body problem. 

The integration for each particle was terminated if it reached 
Too = 200. Integrating to higher values of r„o, up to ~ 3 x 10 , 
yielded negligible effects on the results presented below. An up- 
per limit on the physical distance of the order of < 10 pc is also 
expected as tangential forces produced by various deviations from 
the pure spherical symmetry, assumed here, become more likely 
and more effective in changing the particle's angular momentum, 
when moving far away from the centre. 

The integration was also terminated if the particle reached 
''tidal = lO^"' from either BHs. This represents the tidal disruption 
of the star by the BH (the orbit is than termed as "crashed"). The 
true tidal disruption radius is r»(2/WBH/"'*)''^, where r^, is the ra- 
dius of the star, which is about two orders of magnitude smaller 
than assumed here. However, setting rtijai = 10^ was enough to 
set a negligible rate of stellar depletion compared to the rate of es- 
cape from the system (reaching r^), and thus led to a negligible 
effects on the results. Reducing further rtijai was also expensive in 
computing time, and was therefore avoided. 

The integration was usually stopped at f = lO"*, which cor- 
responds to ~ 800 revolutions for a g = 1 binary, or a physical 
time-scale of ~ 10^ years. Orbits that have neither diverged (es- 
caped) nor crashed before the end of their integration were defined 
as "stable". We also investigated the effects of extending the orbital 
stability to ? = 10 , i.e. extending to 10 years. As we show be- 
low, most unstable orbits diverge on a few orbital time-scales, and 
a steeply decreasing fraction of the phase space volume is popu- 
lated by unstable orbits with increasing time-scale. 
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3.2 Simulations 

To better understand the effects of a BBH on the LOS VDs, we first 
explore the evolution with time of the regions in velocity phase 
space populated by stable orbits, as a function of distance from the 
BBH. For this purpose we use stability maps (see §4.1). In these 
maps, a fixed spatial position is chosen as the particle's starting 
point, and a 2D grid of initial velocity vectors in a certain plane 
is created. Each grid point is integrated up to f = lO'*, or f = 10*, 
and if the integration is terminated earlier, the grid point is tagged 
according to the time it took for the orbit to terminate (reach R„ 
or rtidai). The stability maps are also useful for understanding the 
effects of various parameters (e.g. the bulge potential) on the or- 
bit stability, and are also a tool for the code development (locate 
inaccuracies and bugs throughout parameter space). 

The stability maps do not produce directly observable results. 
For this purpose we performed 3D Monte Carlo (MC) simulations 
to derive the observable kinematic signature. The initial positions 
of the particles in these simulations were drawn randomly from a 
p(r) oc r^~ distribution up to rmax = 60. It was verified that parti- 
cles which were initiated at r,nax > 60 (using a shell 60 < r < 80) 
made a negligible change in the results. The decreasing effect of 
outer regions is expected since the line of sight integration of the 
r stellar distribution scales as r^'. For the sake of simplicity we 
did not introduce a flat core at r < /? = 1 to the initial MC re- 
distribution, as the r < 1 region is only relevant for the spatially 
unresolved stellar population within the Hill sphere of each BH. 

In velocity space, initial distributions in each direction were 
drawn from a normal distribution with variance CJ-, or equivalently 
a Maxwell-Boltzmann f{\'), as noted in §2.3. A position dependent 
cutoff was applied so that |v| < I'esc, where Vesc is the velocity re- 
quired to reach roc from a given point in space (therefore, the initial 
f(v) varies slightly with position). This cut is only for the sake of 
convenience, as all orbits with |v| > Vesc are found to be unstable 
(see the stability maps in §4.1). 

Orbits residing within the Hill spheres of the two BHs, with a 
total energy below the minimal potential energy at the Lagrangian 
points (as measured in the corotating system), cannot escape. All 
these orbits are therefore bound, but stars may be destroyed through 
tidal disruption (see Chen et al. 2009), which is beyond the scope 
of this study (see §3.1). Because of the high accelerations in this 
region, the integration time of a single orbit can exceed by a few or- 
ders of magnitudes the integration time for an orbit outside the Hill 
spheres. We therefore separated each simulation into two parts: one 
for initial conditions within r < 2 which includes almost all the Hill 
sphere orbits, and the other for 2 < r < 60, which includes almost 
none. The two datasets were pieced together with the appropriate 
weights for an r^^ distribution. 
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Figure 1. The evolution of the stability maps with time, up to ( = 10*, for 
a cj = \ BBH. The horizontal and vertical axes coiTespond to the initial v^ 
and Vy velocities, normalized by the local I'esc (where I'esc = 1.31,1.06 for 
X = 5, 10). The orbits are launched from a point along the j:-axis: x = 5 
(top panel) and .r = 10 (bottom), with \\ = (planar orbits), i',, is therefore 
the tangential velocity, and positive values of which belong to corotating 
orbits. The colour indicates the survival time of each orbit. White represents 
stable orbits (survive lo t ^ 10*). The black lines bound the orbits within 
the classical loss cone. Note the near convergence of the stable areas for 
t > lO'. Below we integrate lo t = 10"*, which con'esponds to a ~ 10% 
overestimate of the area in phase space occupied by stable orbits. 



4 RESULTS 

4.1 Stability Maps 

Figure 1 presents cuts in phase space for g' = 1 . Each point (or pixel) 
represents the initial velocity conditions of an orbit. A white pixel 
represents an orbit that neither diverged nor crashed throughout its 
integration; the orbits were followed to f = 10 . Non- white pixel is 
an orbit that became unstable at a time indicated by the colourbar. 
The upper and lower panels represent particles launched from x = 
5, 10, respectively; in all panels {y,z) = (0,0). The cuts in velocity 
space are in the v_^-Vy plane, for orbits with i'- = (i.e. purely planar 
orbits). A total number of 256" orbits were calculated in each map. 



where in each axis grid points are uniformly spaced in the range 
^^'esc < Vi < Vesc- Note that Vgsc differs depending on the launching 



radius (see caption). 

All particles with v^ = v^ + v^ + v^ > v^^^ are unstable, as in- 
dicated by the circular boundary of the stable region, which has 
a radius of unity (Fig. 1). Note that some orbits with v > Vesc, 
which start inwards (vj < 0), do not escape immediately, in con- 
trast to the outgoing orbits (v.v > 0). The orbits which start inwards 
with v > I'esc can be temporarily trapped by the BBH, and wonder 
around on chaotic orbits. However, chaotic orbits are inherently un- 
stable, and they inevitably lead to escape or a tidal disruption on 
time-scales of 10-10 for these orbits. 

The solid black lines correspond to the initial conditions re- 
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quired to reach a pericentre distance , Tperi, of rmjn = 1, assuming a 
purely central force (i.e. the two BHs are taken to be a single point 
source; the bulge force is unchanged). The line shape is given by: 



:± 



2 [*(x) -*(,•„;„)] +v2 



(^/'■ir 



-1 



(17) 



where 'I>(r) = '5bulee('') ^ (1 +?)/'' is the total potential, and the 
bulge term is given by equation (14). 

The orbits arrive to ~ r^in, interact strongly with one of the 
two BHs, and are flung out. For negative values of Vy, correspond- 
ing to retrograde orbits (i.e. opposite to the BHs' rotation), the 
black lines fit well to the boundary of the stability region at f = 10^ 
(especially for the larger r, where the tangential part of the force is 
small). A similar boundary is seen at positive v,., but corresponds to 
higher angular momentum than expected from the simplified loss 
cone solution. 

As the integration time increases, from / = 10 to 10^, the 
boundary of the stable regions in Fig. 1 contracts. At x = 10 there 
is a roughly uniform reduction in the stable area with each ten fold 
increase in time, indicating the reduction in area drops logarith- 
mically with time. At x = 5 there is a similar trend in the area of 
the stable retrograde orbits, but the prograde orbits area is roughly 
stable beyond f = 10^. The simulations described below were car- 
ried out to f = 10 , as the results get reasonably close to conver- 
gence on this time-scale. The asymmetry, which produces the BBH 
kinematic signature discussed in this paper, increases somewhat on 
time-scales longer than 10 . If the system is allowed to evolve fur- 
ther, the kinematic signature will be somewhat enhanced with re- 
spect to the results for f = 10 , presented below. 

In Figure 2 more phase space cuts are presented, for orbits 
followed to f = 10^. The greyscale indicates the time of instability 
while white pixels are stable. The left, middle, and right columns 
represents particles launched from x = 3,5, 10, respectively; in all 
panels {y,z) = (0,0). The upper row shows cuts in velocity space 
in the Vx-Vy plane, for orbits with v- = (purely planar orbits); the 
lower row shows cuts in the Vy-V; plane, for orbits with v.v = (ini- 
tial velocity tangential). Here a blue line indicated the borders of 
the classical loss cone region. 

With decreasing launching radius, the radial orbits get more 
depleted, and the tangential orbits develop a growing asymmetry. 
At X = 3, only retrograde planar orbits remain. 

These plots demonstrate that the "loss cone" term, coined by 
Frank & Rees (1976) is inaccurate. The unstable orbits occupy a 
loss cylinder, as pointed out by Cohn & Kulsrud (1978). The loss 
cylinder grows asymmetrically towards the centre, and dominates 
most of the phase space volume at r < 5. 

The least stable orbits, apart from those on direct collision 
paths with one of the BHs' tidal radii, are the retrograde orbits just 
outside the boundary of the classical loss cone, noted by the lower 
blue lines in Figs. 2 and 3; these orbits appear as dark grey pixels in 
the stability maps, as they escape on short time-scales. These orbits 
lead to head on collisions with the approaching BH, and the stars 
are flung out immediately. However, once the pericentre distance 
of the retrograde orbit is rperi > r^in (orbits below the lower blue 
line), i.e. the orbit is outside the classical loss cone, the BHs cannot 
deflect the star appreciably, and the orbit becomes stable on long 
time-scales. Torquing by the binary on the star changes rapidly due 



^ We use the term pericentre distance to describe the radius of the closest 
point of a star from the centre of the coordinate system. 



to their large relative velocity with the star, and the net torque tends 
to cancel out in each close approach of the star. 

In contrast, prograde orbits close to the classical loss cone are 
flung out on longer time-scales. In this case the star approaches the 
receding BH and the relative velocity is small; the star is subject to 
a nearly steady torque by the binary, which leads to some energy 
gain. Repeated close encounters build up the energy of the star, 
until it is ejected. The further away the pericentre distance is, the 
smaller is the energy gain, and the longer it takes the star to build up 
the ejection energy. This leads to a gradual increase in the escape 
time, moving away from the prograde classical loss cone boundary. 
Another possible scenario is that prograde orbits close to the clas- 
sical loss cone do not gradually build up their energy, but rather are 
on quasi-regular orbits with a nearly fixed energy, i.e. on chaotic 
orbits with a longer Lyapunov time, which increases with distance 
from the centre. Once the orbit becomes chaotic, a disruptive head 
on collision with one of the BHs is quick and inevitable. 

The same concept is also illustrated in Figure 3, but there 
5 = 0.1 and thus r^^ = i?2 ~ 1-8, the orbital radius of the secondary 
BH (see equation 5). The larger r^{^ in the ^ = 0. 1 case, leads to a 
larger loss cone region, relative to Vgsc, and thus to a smaller area in 
phase space available for stable orbits. Thus, a g < 1 binary leads 
to a larger effect on the orbits stability compared to a <? = 1 binary. 
However, in the limit g — >■ the effect of the secondary clearly dis- 
appears. 



4.2 The Projected Kinematic Signature 

4.2.1 Monte Carlo Simulations 

Four MC simulations were carried out to calculate the projected 
LOSVDs. To increase the statistics for the kinematical maps shown 
below, we took a number of snapshots of the system at different 
times; all snapshots were taken at times in which the BHs had the 
same orbital phase (position on the x-axis). The snapshots were 
taken close to the end of the integration at f = 10"*, where the sys- 
tem was closest to a steady state solution, and were spaced in time 
by one BBH revolution to ensure significant offsets of the stars be- 
tween snapshots. 

The simulations main and main2 are for g = I. In main, the 
binary is observed along the y-axis (edge-on) and the BHs are max- 
imally separated to the observer (side-view). In main2 the view 
is along the x-axis so the BHs are along the line of sight (front- 
view). In these two simulations a = 0.25 (see Table 1). Simulation 
ratio is for g = 0.1, for an edge-on side-view. In this calculation 
(7 ~ 0. 1 1. The simulation single is for a single BH, located at the 
origin (i.e. q = 0); we used CT = 0.25 for it. Note that for this case, 
the binary separation a has no meaning, and thus there is another 
free parameter to determine the physical scales of the system. 

Table 2 lists for each simulation the total number of orbits 
integrated, where A' (inner) and N (outer) refer to orbits initiated 
at r < 2, and 2 < r < 60, respectively (see §3.2). The values of A' 
were set to be large enough to minimize the statistical noise in the 
results, and may be well above the number of stars contributing to 
the observed stellar spectral features in typical galaxies. 

Table 2 also gives the fraction of the orbits which are classi- 
fied as divergent (reaching r„) or crashing (reaching rjidai). This 
fraction was calculated by giving the weights to the inner and outer 
datasets, expected from the r^ density profile. The q = 0.1 simu- 
lation has ~ 30% more diverging orbits, and twice the fraction of 
crashing orbits than the q = I simulation. However, the total bulge 
masses of the models are different (depending on a), and there- 
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Figure 2. Stability maps for cuts along different planes, as a function of distance from the centre. In both rows, the orbits are started on the x-axis and are 
launched from .r = 3 (left column), ,v = 5 (middle) and x = 10 (right). The integration time is ? = 10**, for aq= \ BBH. The upper panels represent orbits in 
the x-y plane (iv = 0), and the lower panels orbits start in the y-z plane {v., = 0), i.e. start as purely tangential. Greyscale indicates the survival time of each 
orbit, and the white areas are stable orbits The panels show the initial velocities normaUzed by Ves^ (= 1-57, 1.31, 1.06 for x = 3,5, 10). Note the decreasing 
area of stable prograde orbits (v,, > 0) with decreasing distance from the binary. Almost only retrograde orbits are present at x = 3. The blue Hues represent 
the border of the classical loss cone. It reproduces fairly well the boundary of the stable retrograde orbits, in particular in the x-y plane, but fails significantly 
for the prograde orbits, already at ,v = 10. The dashed red hne represent the fixed a of the isothermal Maxwell-Boltzmann f{v) used below. 



fore these fractions translate into different mass deficits. For the 
q = \ simulations, the total bulge mass at rmax = 60 was M ~ 7.42 
(see equation 13). The relative mass deficit due to diverging orbits 
in these simulations is therefore Mi^f/Mi2 ~ 1.01, where M12 is 
the combined BH mass, l-\-q. For the ^ = 0.1 simulation, the to- 
tal mass at rmax was M ~ 1.35, and the mass of diverging orbits 
was M(ief/Mi2 ~ 0.45. The larger fraction of diverging orbits for 
g = 0. 1 results from the lower bulge contribution within the binary 
orbit, as expressed by the lower a jV , which allowed more orbits 
to reach r > r„. As, expected, the two q= \ simulations yield the 
same fraction of divergent and crashing orbits, as the difference is 
only in perspective, and the statistical error is small enough. 

The fraction of 10^ orbits which diverged in the ^ = sim- 
ulation results from the numerical error in energy conservation, 
which allowed this fraction of orbits with v < I'esc to become un- 
bound. The fraction of crashing orbits is a factor of 3 to 6 times 



smaller than in the ^ > simulations. This demonstrates qualita- 
tively the enhancement of a BBH on the tidal disruption rate (e.g. 
Chen et al. 2008; 2009), though the exact numbers are not valid 
given the high value of rtijai used here. In the ^ = 0.1 simula- 
tion, the primary tidally disrupted 67 times more stars than the sec- 
ondary. 



4.2.2 The Velocity Distribution Moments 

To characterize the shape of a LOSVD using a small num- 
ber of parameters, we expand it to a series of the so called 
Gauss-Hermite (GH) moments; this procedure is consistent with 
van der Marel & Franx (1993). We use a function of the form: 



^(v) 



y-i^ 



271(7 



N 



(18) 



© 2010 RAS, MNRAS 000, 1-18 



The Stellar Kinematic Signature of Massive Black Hole Binaries 



-1.0 -0.5 0.0 0.5 



-1.0 -0.5 0.0 0.5 



1.0 -0.5 0.0 0.5 




10- 



10-' 



10" 



Figure 3. Same as Fig. 2 for a ^ = 0.1 BBH. Axes are normalized by v'esc (= 0.96,0.77,0.59 for x = 3,5, 10). Note that despite the factor of 10 decrease in the 
secondary BH mass compared to the q=\ case, the loss cone in fact gets larger (see text), and significant progradre/retrograde asynmietry remains. 



where w = (v — /i)/c7 is the normalized velocity parameter, H„{w) 
is the Hermite polynomial of the «* degree. Note that CT here refers 
to the dispersion of this particular line profile. We find the best- 
fitting y, /x, a and h„ {n ^ 3) parameters using the least-squares 
method. 

While the GH moments derived above are not standard mo- 
ments in the statistical sense, as they are derived by least-squares 
best-fitting to the data, rather than by projections on the data, these 
commonly used GH terms are useful to describe the deviation of 
the velocity profile from a Gaussian. For practical reasons we are 
mostly interested in the first two deviations, as these are often well 
measured. The corresponding polynomials are: 



«3W 



1 

7! 



Ix" - 3x 



HAx) = — - (Ax* - 12x- + 3 
2^/6 V 



(19) 
(20) 



The coefficient of H^{x), Ag, is a measure of the profile's lack of 
symmetry for reflections with respect to its centroid, and /14 is a 
measure of the even deviation from the Gaussian shape. However, 
there can be significant contributions to the profiles from higher 



moment and we therefore also present below plots of the actual 
LOSVDs along various positions. 



4.2.3 Projected views 

Figure 4 presents maps of the projected ^, a, /13 and h^ for the q=\ 
simulation; for an edge-on side-view (along the ^^-axis), where the 
BHs are maximally separated. Each map is composed of 49 x 49 
pixels in the projected plane; each pixel represents a moment calcu- 
lated from the LOSVD measured using 100 bins in velocity. Thus, 
the spatial resolution is ~ 0.4 length units. Due to the maps' up- 
down symmetry and left-right antisymmetry, and the the fact that 
each is made of 10 superimposed snapshots, statistics is increased 
40-fold. Therefore, these maps are derived from a total effective 
number of ~ 3 x 10 particles. 

A prominent counter rotating "torus"-like structure is seen in 
the top-left /x map out to a scale 5 times larger than the binary 
separation. This results from the preferential stability of retrograde 
orbits, clearly seen in the stability maps (Fig. 2). On scales of r < 2 
one can see the prograde orbits of the stars trapped in the BHs' Hill 
spheres, as they move together with the BHs. 
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Figure 4. Maps of velocity moments of the LOSVD for a ty = 1 BBH. An edge-on side-view of the BBH. The BHs are located at x = ±\,{y,z) = (0,0). 
The line of sight velocity corresponds to v,.. Note the "torus"-like structure in n produced by the dominance of retrograde orbits outside the BH Hill spheres. 
The prograde motion within the Hill spheres reflects the initial conditions. There is a dip in CT towards the centre, in contrast to the monotonic rise around a 
single BH. A torus and a dip are present also in the maps of the higher moments A3 and /14. Note that the peak values of the pixels within the Hill spheres ai'e 
saturated. 



The top-right panel in Fig. 4 shows a map of a. As expected, 
a increases inwards, and peaks inside the Hill spheres of the two 
BHs. However, there is a prominent drop in CJ at r < 3. This drop 
occurs at r where all the stable orbits are purely retrograde, and are 
thus moving relatively coherently around the BBH. Inside the Hill 
spheres, both prograde and retrograde orbits are allowed, and the a 
jumps to the expected value for the kinematics around a single BH. 
The maps of /13 and h^, on the lower left and lower right panels, 
show similar structures to those seen in the /J and o maps. 

Figure 5 shows the projected kinematics for an edge-on, front- 
view (BBH viewed along the binary axis). The maps are remark- 
ably similar to the side-view maps on scales r > 2. On smaller 
scales only a single peak is seen in CT, as the two BHs project on the 
same position. Interestingly, the /x map still shows a prograde struc- 
ture, although the two BHs are moving tangentially to the line of 
sight. This may result from the non-uniform and anisotropic initial 
conditions of the bound orbits within the Hill spheres, as measured 
in the local corotating frame centred on each BH. 

Figure 6 shows a face-on projection of the kinematics (the bi- 
nary is viewed along the z-axis). From symmetry, /x = and /i3 = 



everywhere, and only the maps of a and h^ are shown. Again, a 
similar structure is seen on scales r > 2 to that seen in the above 
two perspectives. The small difference is the perfect axial symme- 
try, in contrast to the reflection symmetry in the edge-on views. 

Figure 7 presents an edge-on side-view of kinematics around 
the ^ = 0.1 BBH. The more massive BH is on the right. The 
kinematic signatures remain prominent, but the structure obviously 
loses the reflection symmetry of the q= i case. Also, the maps of 
/73 and /74 become more distinct compared to the maps of jl and a, 
and not similar as they were in the q = \ case. The face-on view 
in Figure 8 shows that the low mass companion creates a low CT 
"trench" along its orbit, and a tiny peak within its tiny Hill sphere. 
A similar effect is seen in the h^ map. 



4.2.4 Slit Views 

Figure 9 presents slit views of the four velocity distribution mo- 
ments. The slit is placed along the x-axis; and provides a cross sec- 
tion of the edge-on side-view kinematical maps. The solid lines is 
for the q = \ case, and the dashed line for the q = Q.\ case. The 
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Figure 5. The same as Fig. 4, for an edge-on front-view, i.e. along the ,v-axis, so both BHs are along the line of sight. As expected, the structure on scales 
beyond r>2 remains the same. The Hill spheres now overlap and thus form a more compact structure. 
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Figure 6. The same as Fig. 4, for a face-on view (along the ;-axis). Due to the reflection symmetry of the system with respect to the x-y plane, the LOSVD 
is symmetric, and thus fl = h^ = 0. The structure in the maps of c and /14 outside the Hill spheres is similar to that seen in the edge-on views, but now the 
structure shows a perfect axial symmetiy, in contrast to the reflection symmetry in the edge-on views. 
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Figure 7. An edge-on side-view of aq = 0. 1 BBH, as in Fig. 4 for ij = 1. The primary BH is now at j = 0.18 and the secondary atx = — 1.82. The structures 
seen in the ll and c maps are similar to those seen in the q=\ case, but a is now physically 5.5 smaller for the same M, (Table 1). The amplitude of jl in the 
torus structure, is now lower, in model units, compared to the q=\ case, but this is compensated by the higher physical value of the velocity unit. The higher 
GH moments display more complicated structure, but roughly similar to the ^ = 1 case. 
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Figure 8. A face-on view of a </ = 0. 1 BBH, as in Fig. 6 for aq=\. The secondary BH at .t = —1.82 now scours a "trench" in fj and /!4, rather than a wide 
dip seen in the q= \ case. 
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Figure 9. The view through a slit of the edge-on side-view maps presented in Figs. 4, 7. The slit is placed along the x-axis and extends from z = —0.6 to 
z = 0.6. The g = result is scaled to M, = 2, so a comparison with the q=\ can be made. Since /(v) remains isotropic for the i^ = case, both \i. and /13 
are there. Note the higher CT at |.v|> 3 for the q=\ case, compared to a single BH of the same total mass. Also, the drop in CT at |.r|< 5, in contrast to the 
q = Q case. Lai'ge amplitude features are also seen in /13 and hi,, in contrast to the g = case. Note that the horizontal and vertical axes correspond to different 
physical scales for different values of q, for the same M, (see Table 1 ). 
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Figure 10. The same as Fig. 9, but the system is viewed here along the z-axis ("face-on"). The odd moments, ;U and /13, are not shown as they are zero due to 
symmetry. 
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Table 2. The properties of the four MC simulation carried out. N is the 
number of orbits integrated in each run (outer and inner regions were calcu- 
lated separately, see §3.2). The number of snapshots used in the Figures of 
§4.2 is indicated in the fourth row. The two bottom rows give the fraction 
of divergent (r > r^) and crashing (r < rudai) orbits. These numbers can be 
converted to relative mass deficiencies (see §4.2.1). 
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top left panel present /x. The two sharp peaks at .v = ±1 for q= 1 
represent the cluster of stars trapped in the Hill sphere of each BH, 
moving at approximately the BHs orbital velocities of ±0.5. Fur- 
ther out, at X ~ ±3, there are the broad peaks of the larger scale 
counter rotating torus, as seen in Fig. 4. At a low angular resolu- 
tion, the broad counter rotating torus structure may be masked by 
the compact corotating clusters. The amount of dilution depends 
on the compact clusters luminosity compared to the torus stellar 
luminosity, which may be different than assumed here. 

Interestingly, when q = Q.i, the two sharp peaks disappear. 
The massive component, located &X x = R\ ~ 0.2, now moves at 
a velocity of only I'l = qj ^2{\ +^) r^ 0.067, below the torus 
peak velocity. The secondary BH moves faster now, with i'2 = 
1/1^2(1 +(^) ~ 0.67, but its Hill sphere volume is now roughly 



7-3/2 ^ 



: 30 times smaller, and its contribution to the total profile. 



given the assumed p(r), is negligible. 

The top right panel of Fig. 9 presents a along the slit. The 
dotted line is for a </ = single BH, with a total mass scaled to 
2, which shows the expected monotonic rise as a oc ^-1/2 towards 
the centre. In the q=\ case, C7 at x > 3 is larger by 20-40% than 
for the single BH case, with the same total mass. This rise occurs 
because of the exclusions of the low velocities within the loss cone. 
Furthermore, at x < 5, a starts falling, in contrast to the sharp rise 
for the ^ = case. This results from the gradual elimination of 
the prograde orbits with decreasing r. This leaves only retrograde 
orbits, and a quasi-coherent flow, and thus a lower a. The double 
peaks at the centre are due to stars bound within each Hill sphere. 
For g = 0. 1 only one peak is prominently seen, as expected, since 
the volume of the secondary Hill sphere is ~ 30 times smaller. The 
radial profile of CT also shows a small excess compared to the single 
BH with the same total mass"*, and a drop close to the centre, but 
the effects are less pronounced then for the q=\ case. 

The lower left panel shows /13 along the slit. Comparable peak 
values are seen in both the q=\ and ^ = 0.1 cases, indicating that 
the amplitude of the LOSVD asymmetry is driven by the presence 
of a second BH, and is not sensitive to its mass for q = m the range 
0.1-1. This can also be seen in the stability maps, which show sim- 
ilar asymmetry for q= \ and g = 0.1. The value of h^ along the 

** the plotted q = can be scaled to a total mass of 1 . 1 , for the q = 0.1 case, 
by multiplying the velocities by \/2/ll a: 0.43, and distances by 121/40 ~ 
3. 



slit is shown in the lower right panel. The differences from the sin- 
gle BH case are more prominent, and q = 0.1 presents the largest 
effects. 

Figure 10 presents the slit results for a face-on view. The ex- 
cess in (7, and its depression at \x\ < 4 is clearly apparent for the 
q = i binary, and a somewhat enhanced depression is also seen for 
the q = 0.l binary, compared with the edge-on view. Significant de- 
viation of the h4 profile from the single BH are now prominent for 
both q values. 



4.3 Density Profile 

Figures 1 1 and 12 present an edge-on side-view, and a face-on view 
of the surface stellar densities for the q= I and ^ = 0.1 cases. Both 
images for the q = I case are remarkably similar, showing the well 
known "scouring" effect of the BBH which depletes stars close to 
the binary. The slight difference is the perfect axial symmetry of 
the core structure in the face-on view versus the somewhat elon- 
gated core structure in vertical direction in the edge-on view. In 
the ^ = 0.1 case the edge-on side-view shows that the secondary 
BH carves out a low density torus structure along its orbit around 
the primary, producing a circular structure extended in the vertical 
direction, for the innermost stellar light distribution around the pri- 
mary BH. The face-on view shows an axially symmetric structure 
with a tiny density enhancement from stars bound to the Hill sphere 
of the secondary BH. 

Figure 13 shows the radial density profile p(r) averaged along 
spherical shells. The MC initial condition is p{r) oc r-" out to 
r = Rmax = 60, as noted by the red line in the figure. In the ^ = 
simulation the slope remains close to —2 in the inner parts. The 
slight steepening of the slope towards Rmax, and the much steeper 
slope beyond Rmax, are edge effects which come from the fact that 
particles on radial orbits from the initial r < 60 sphere move out to 
60 < r < roo. 

For the BBH cases, the slope flattens inwards towards the core, 
as the fraction of unstable orbits increases with decreasing r. In the 
q = I case, a dip forms for r < 5, with a minimum at r ~ 2. The 
rise inward at r < 2 is produced by stars bound within the Hill 
spheres of the two BHs. In the q = 0.1 case the dip is shallower, as 
the averaging over spherical shells dilutes the density drop which 
occurs only in the region close to the secondary BH. The difference 
between the initial and final density distributions is a measure of the 
stars lost from the system as a function of r, as noted in §4.2.1. It 
is important to note that although our initial condition has a core 
radius h= I, the final state core radius is ~ 10, as shown in Fig. 13, 
which is closer to the observed core radii in core ellipticals. 



4.4 Internal Kinematics 

Figure 14, upper panel, shows the mean tangential component of 
the velocity, V,;,, as a function of r. At r > 10 both prograde and 
retrograde orbits are similarly stable, and thus V^ approaches zero. 
At r < 10 retrograde orbits become significantly more stable than 
prograde (Fig. 2), leading to a sharp rise in V^. For q = l,V^ rises 
from ~ 10% of the circular velocity at r = 8, to ~ 50% at r = 4 
to ~ 100% at r = 2, where only retrograde orbits survive. The 
solid red line is the circular velocity, for comparison; it was cal- 
culated assuming a BH mass of 2 in the centre, as in the q = 1 
case. For q = 0.1, the asymmetry between retro- and prograde or- 
bits remains significant, as can be seen in Fig. 3, but averaging over 
spherical shells dilutes the effect of the secondary BH, due to the 
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Figure 11. The projected density for the q = I simulation. The left panel is an edge-on side-view, and the right panel is a face-on view. Both images are 
remarkably similar, showing the well known "scouring" effect related to the BBH formation. The slight differences is the perfect axial symmetry in the face-on 
view, and the somewhat elongated core structure in vertical dii'ection in the edge-on view. Colour represents the logarithm of the normalized density. The 
values in the inner regions ai'e saturated. 
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Figure 12. The same as Fig. 11, for a ^ = 0.1 BBH. The secondaiy BH carves out a low density torus structure along its orbit around the primaiy, producing 
a circular structure extended in the vertical direction for the innermost stellar light distribution around the primary BH. The face-on view shows an axially 
symmetric structure with a tiny density enhancement from stars bound to the Hill sphere of the secondary BH. 



non-axisymmetric morphology of the velocity field in this case, as 
evident in Fig. 7. 

The middle panel of Fig. 14 shows the 3D velocity dispersion 
ct' = [(cT^^ + O'i + CT^)/3] '/^ (not to be confused with the ID veloc- 
ity dispersion, a). The flattening and drop at r < 3 results, as noted 
above, from the disappearance of prograde orbits, which leads to 
a more coherent flow with only retrograde orbits, and thus a lower 
dispersion. The red line indicates a' for the ^ = simulation, for a 
single BH with a mass of 2. The excess of the BBH ID ct, noted in 
the slit view, can also be seen here for the 3D ct'. 

The bottom panel of Fig. 14 shows the anisotropy parameter, 
J3 = 1 — a} /o}, where ct, = (ct? + CTg)''~ is the tangential veloc- 
ity dispersion, and CTr is the radial velocity dispersion. The velocity 
field is significantly anisotropic already at i?,nax> as orbits within the 
loss cone are excluded, leading to ct, > a,- throughout the shown r- 
range. As the loss cone grows inwards, the velocities become more 



tangential, and thus /3 becomes more negative. At r < 3 only retro- 
grade orbits remain and the tangential orbits become more coher- 
ent, which reduces ct, and increases /3 inwards. 

The anisotropy derived here is about an order of magnitude 
larger in absolute value than in Milosavljevic & Merritt (2001) (see 
their figure 16). We suspect that this issue results from the fact that 
the BBH in that work does not stall at «[, but rather still shrinks 
rapidly at this radius (as indicated by figure 1 there). The stellar dis- 
tribution there cannot reach a near steady state solution, as derived 
here. We verified that a shorter integration time in our simulation 
indeed yields a j3 which is smaller by factor of few. Thus, the stel- 
lar kinematics derived in Milosavljevic & Merritt (2001) does not 
capture the full level of the kinematic signature of a BBH which is 
stalled at aji- 

The dotted line model, shown in all panels, presents a modified 
q= i model with a uniform velocity distribution extending to Vgsc, 
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Figure 13. The radial stellar density distributions for different q values. 
n(r) is the normalized number of particles per unit volume. The red line is 
the initial « r^^ distribution in the r = 1-60 range. In the q = case the 
distribution remains nearly unchanged, apart from an extension to r > 60 
by stars on nearly radial orbits. The q = I model shows the flattening of 
p{r) at r < 5, and the deep minimum at r = 2, just outside the Hill spheres 
of both BHs. In the q = 0.l case the dip is significantly reduced, partly due 
to the averaging over spherical shells of a stronger dip confined only close 
to the secondary BH orbit (see Fig. 12). 



or effectively a Maxwell-Boltzmann distribution with CT — > oo. The 
remarkable similarity to the standard q = I model, where a = 0.25 
(Table 1), for the three internal kinematics parameters, shows that 
the BBH kinematic signature is nearly independent of the form of 
/(v) used for the initial conditions. 



4.5 LOSVDs 

Figure 15 shows the LOSVD along certain lines of sight. Each 
panel shows the LOSVD along three lines of sight at distances of 
3, 5, and 10 from the centre of gravity (origin of the coordinates). 
The size of the aperture corresponds to a projection pixel, an area 
of 0.4 X 0.4 square length units. The top panels are for the q = I 
case, and the bottom panels are for q = 0.1. The left panels are for 
an edge-on side-view of the binary and the right panels are for a 
face-on view. The vertical axis is in units of orbits per velocity bin, 
where each line has 100 velocity bins; the range of velocities is de- 
termined according to the maximal escape velocity along each line 
of sight, the same values used in the stability maps (Figs. 2 and 3). 
The line FWHM changes by ~ 10% moving inward from 
.z = 10 to X = 5, for the q = I model, and it drops moving further 
inwards to x = 3, as also seen in the slit views of a (Figs. 9, 10). 
This is in sharp contrast to the single BH case, where the FWHM 
is expected to rise by 83% (= ylO/S) from x = 10 to x = 3. Note 
also the line asymmetry in the edge-on view, which increases mov- 
ing inwards, reflecting the enhanced retrograde motion close to the 
centre. In the ^ = 0.1 case the presented lines of sight are away 
from the secondary BH, at x = — 1 .82, somewhat reducing the pro- 
file asymmetry. The FWHM of the lines increases from x = 10 to 
X = 5, as expected for a single BH, but it remains constant from 
X = 5 to X = 3 (as seen in Figs. 9, 10), in contrast to the expected 
30% rise for a single BH. A noticeable asymmetry near the line 
base is seen for the edge-on view, similar to the asymmetry seen in 
the q = I case, but with a lower amplitude. 



00. 




Figure 14. The radial dependence of the three internal kinematics parame- 
ters. The top panel is the mean tangential velocity. At r < 10 the preference 
for retrograde orbits becomes significant, leading to the rise in V.^,, reach- 
ing a pure circular velocity (iv) at r < 3, as only retrograde orbits remain. 
The middle panel presents the 3D rms velocity dispersion. At large r the 
q= I model matches the ^ = case (where M, = 2), but at smaller r there 
is a 20-30% excess (see also Fig. 9) produced by the BBH modification 
of f{v). The bottom panel is the anisotropy parameter /3. The q = I and 
q = O.l models show the velocities become tangential close to the cen- 
tre, in shaip contrast to q = 0, where /(v) remains nearly isotropic at all 
r. The model u.v., shown in all panels, presents a modified q = 1 model 
with a uniform velocity distribution extending to I'esc (effectively Maxwell- 
Boltzmann with (T — > oo). The remarkable similarity to the standard q = 1 
model where a = 0.25 (Table 1), for all three parameters, shows that the 
BBH kinematic signature is nearly independent of the form of f{v) used for 
the initial conditions. 



The LOSVDs for V;, seen for the face-on view, are symmetric, 
as expected due to the reflection symmetry of the system with re- 
spect to the x-y plane. In contrast, the edge-on v,. profiles are asym- 
metric due to the prograde/retrograde asymmetry discussed above. 
The asymmetry increases for lines of sights closer to the centre, 
and towards the wings in each profile, as these are produced by or- 
bits closer to the binary, where the tangential velocity asymmetry 
becomes larger. 

Figure 16 shows the LOSVDs expected from a low angular 
resolution observation for an edge-on line of sight. We compare the 
velocity profiles for two lines of sight situated on opposite sides of 
the centre, at distances of 5 and 10 from the centre. We integrated 
the LOSVD through a circular aperture using a Gaussian with a 
FWHM of 10 as the weight function, which represents the angular 
point spread function (PSF) of the telescope. The left panel is for 
q = I and the right for g = 0.1. The q = 1 case shows a clear shift 
in the peaks of the LOSVDs, which results from the net rotation of 
the stars, clearly seen in the spatially resolved maps of the average 
velocity fl (Figs. 4, 7). The shift is more prominent for the lines 
of sights centred at x = ±10, i.e. at a position a FWHM of the 
PSF away from the centre. In the ^ = 0.1 case the PSF used here 
eliminates almost completely the profile differences between the 
two lines of sight, and the binary kinematic signature is very small. 
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5 DISCUSSION 

The extensive scattering experiments presented above, using ~ 10 
test particles surrounding a massive BBH embedded in a bulge po- 
tential, allowed us to accurately map the 3D velocity distribution 
of stable orbits. These are used to derive maps of the projected 2D 
velocity distribution moments, and of the LOSVD along various 
directions. The stable orbits close to the binary are generally tan- 
gential and preferentially retrograde, leading to retrograde "torus" 
structure in the projected average velocity. The velocity dispersion 
shows an excess of 20-A0% compared to a single BH of the same 
total mass, and shows a dip close to the binary. These effects lead to 
a clear kinematic signature of the BBH, which can be detected on 
scales of 5-lOa. Thus, they can be spatially resolved even when the 
binary cannot be resolved, as expected even in the nearest galaxies 
(Yu 2002). 

Interestingly, the maps of the 2D velocity distribution mo- 
ments for both the ^ = 1 and the q = 0.1 cases, show that the 
kinematic signature extends on similar scales (in units of r/a) and 
with similar amplitudes (in absolute velocity). One could expect a 
smaller effect if the companion BH is of smaller mass, but appar- 
ently the stability of orbits is strongly influenced by a secondary 
with just 10% of the primary's mass. However, the stalling radius 
is a factor of (1 +q)/2q = 5.5 smaller in the latter case, and will 
thus be harder to spatially resolve. 

The tendency for counter rotating orbits and a veloc- 
ity dispersion drop close to the centre were briefly noted by 
Milosavljevic & Merritt (2001). However, their results were based 
on W-body simulations of ~ 10^ particles, carried out on scales 
'-^100 larger than here, required to simulate the merger of the two 
bulges. As a result, there were only ~ 10 particles in their study on 
the rjnfl scale (see their table 2). The implied large statistical errors 
in that study, therefore did not allow to probe the stellar kinematics 
on the scale of r[„f[ and closer to the BBH probed here. The maps 
of the projected kinematics produced by Milosavljevic & Merritt 
(2001) therefore do not show the BBH signature presented here. 

We also find a clear drop in the projected stellar surface den- 
sity, as stars are efficiently ejected from regions just outside the 
Hill spheres of both BHs (for q = 1), or of the secondary BH (for 
^ = 0.1). This is a well known effect (e.g. Zier & Biermann 2001 ), 
and therefore we do not discuss it further here. 

The advantage of an A'-body simulation is that it allows 
to follow the system, starting with plausible initial conditions 
of separate bulges, and derive the resulting stellar velocity dis- 
tributions on large scales following the merger. For example, 
Milosavljevic & Merritt (2001) found that stars in the merged bulge 
have a net rotation in parallel to the initial angular momentum 
of the two bulges. In the scattering experiments, we just assumed 
some initial conditions for the stellar distribution, and do not know 
whether these were plausible. However, since our simulations are 
much faster than A'-body, we can explore the dependence of the 
results on the initial conditions. The initial stellar f{v) used in the 
MC simulation is drawn from the bulge CT, and as shown in Figs. 2, 
and in particular 3, o can become a small fraction of Vgsc at small 
r. To explore the dependence of the BBH signature on the initial 
a chosen for /(v), we repeated the analysis with the extreme as- 
sumption (J — > CO, i.e. a uniform f(v) extending to Vesc- As shown 
in Fig. 14, this uniform initial f{v) produced nearly identical re- 
sults. The uniform /(v) also led to small deviations (< 10%) in the 
slit view results, indicating that the BBH signature is independent 
of the form of the initial /(v). Thus, the merger history of the bi- 
nary, which may set f{v) on larger scales of the bulge, is likely not 



important on the scales of r/a ~ a few, when the binary becomes 
hard and stalls. 

We also explored the effects of the uniform density core ra- 
dius h (equation 12), and produced stability maps for h = 3 instead 
of h = 1, for the same O. The maps look nearly indistinguishable 
from those shown in §4. 1 . This is expected as increasing h reduces 
the bulge mass, which is already very small for h = I close to the 
BBH (equation 13). So, the BBH signature is independent of the 
exact form of the inner bulge potential, as expected since a hard 
binary resides at a;, <C GM,/4(y, where the bulge mass is neg- 
ligible. However, in a flat core the integrated line of sight stellar 
light increases outwards, increasing the dilution of the BBH kine- 
matic signature by the extended stellar light. The final state core 
radius derived here of H ^ 10 (Fig. 13) is significantly larger than 
the initial condition of h = 1, but is clearly too small in some core 
ellipticals, where the BBH kinematic signature will be harder to 
detect. 

Furthermore, as shown in Fig. 1, a large fraction of the unsta- 
ble orbits are lost on a time-scale oft< 10 , or < 80 BBH periods, 
which corresponds to a few orbital times for orbits starting at r ~ 
5-10. Thus, the BBH kinematic signature is largely imprinted al- 
ready on a time-scale < 10 yr. Unstable orbits present in the initial 
/(v) are quickly excluded. The uniqueness of the BBH signature 
depends on the population of the velocity phase space of stable 
orbits. Since the stable orbits around a BBH are quasi-periodic, a 
given orbit likely moves around in velocity phase space, and so a 
population of stars which populate only a small peculiar corner of 
the stable region in velocity phase space appears unlikely. The BBH 
signature is therefore likely well defined. 

The unstable orbits are essentially orbits within the loss cone 
(or "loss cylinder", as mentioned above), which acquires a more 
complicated shape when r/a is a few. The kinematics we described 
may be more accurately termed "kinematic signature of loss cone 
depletion". The loss cone refilling mechanisms will tend to erase 
the signature presented above. These mechanisms tend to become 
more effective on larger r, in particular at r 2> r[^fi, where either 
steady state, or time dependent perturbations to the spherically 
symmetric potential assumed here are more likely to be found, or to 
occur transiently. On the scale of r < 10 probed here, such mecha- 
nisms are less likely to occur If they occur transiently, they are less 
likely to have a significant effect, given the shorter survival time of 
unstable orbits on these scales. 

The enhancement of the observed o within the BBH rjnfi by 
20-40% implies that the standard direct estimate of M,, which as- 
sumes an isotropic /(v), will lead to an overestimate of M, by a 
factor of 1.5-2. 

Once the BHs coalesced, the signature will be lost on a time- 
scale of Mf/M^, where M» is the loss cone refilling rate and M, is 
the stellar mass at r < 10. This will likely occur on a time-scale sig- 
nificantly longer than the dynamical time-scale at nnj, and may in 
fact be longer than the Hubble time, if it occurs in the low density 
core of a giant elliptical galaxy. However, if the BBH acquires a sig- 
nificant kick following the merger, it may oscillate around the core 
(Gualandris & Merritt 2008), which will erase the BBH kinematic 
signature on a much shorter time-scale, possibly while enhancing 
the core due to heating of the stars. The detection of the predicted 
BBH kinematic signature implies the presence of a BBH currently, 
or a merger which took place on time-scales shorter than the loss 
cone refilling time, although this timescale may be the Hubble time 
in giant ellipticals. 

The calculations presented above follow test particles, and 
thus do not take into account the energy and angular momentum 
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Figure 15. The LOSVDs at different positions. The viewing angle and the model are noted in each panel. Each panel shows three lines of sight (along the 
y-axis) at z = and jr = 3,5, 10. The error bars are statistical errors. The solid lines are the GH best-fits up to order 20. Note that the line FWHM changes by 
~ 10% moving inward from x= 10 to ,r = 5, for the q= \ case, and drops moving inwai'ds to x = 3, in sharp contrast to expected rise for the single BH case. 
In the q = 0.\ case the FWHM increases from ,v = 10 to x = 5, but remains constant from .v = 5 to x = 3, in contrast to the expected rise for a single BH. A 
noticeable asymmetry near the line base is seen for both q values. The fines are symmetric in the face-on views, and show similar trends in the FWHM vs. 
distance from the centre, as seen in the edge-on, side-view. 



lost from the BBH due to the stellar ejection, which can be signifi- 
cant given the large fraction of ejected stellar mass (see §4.2.1). The 
justification is that the purpose of this work is to look for the steady 
state solution, i.e. find which orbits may be populated and which 
cannot survive for long, when the BBH is at the stalling radius, 
rather than follow the time evolution of the system. The implied 
significant energy loss of the BBH found here, results from the in- 
appropriate initial conditions of a spherically symmetric p oc r^- 
assumed here for a hard BBH. In reality, stars may be ejected from 
the system much earlier when the binary just becomes bound, i.e. 
when a ~ r[nfl, or potentially even earlier and on larger scales, based 
on the high value of M^^f/M, ~ 10, and the small scatter, found by 
Kormendy & Bender (2009) for massive ellipticals. 

Here we find that the kinematic signature of the BBH is im- 
printed on the same scale (r < 10) that the surface density signature 
of the BBH is imprinted (inevitable as stars with specific kinemat- 
ics are lost). Therefore, if Mjgf is indeed a signature of the BBH 
formation process, then the BBH kinematic signature may be im- 
printed already on the significantly larger scales of the core radius, 



where Mjgf is measured, of the order of tens to hundreds of parsecs 
(Faber et al. 1997; Kormendy et al. 2009), and may be more easily 
detectable, possibly already in existing data. Clearly, it is interest- 
ing to explore the BBH merger starting from rjng, and follow the 
resulting kinematic signature on larger scales than those calculated 
in this study. 

The calculations presented here assume circular BH orbits. 
The BBH may have a high eccentricity due to various mechanism 
(e.g. Makinoetal. 1993; Mayer et al. 2007; Sesanaetal. 2007; 
Berentzen et al. 2009). In that case the effects calculated here will 
likely extend to larger scales, set by the major axis of the binary 
orbit. It is less clear if the binary stalls in such a case, and at what 
radius, if it does. 

It is also interesting to note that stars bound within the Hill 
spheres preserve the original populations before the merger oc- 
curred, as stars outside the Hill sphere with a total energy below 
the Hill sphere potential barrier (measured in the corotating frame, 
where the energy of each orbit is conserved) cannot enter it, and 
stars with a total energy above the potential barrier in the Hill 



© 2010 RAS, MNRAS 000, 1-lS 



The Stellar Kinematic Signature of Massive Black Hole Binaries 17 




Figure 16. The LOSVDs for an edge-on view with a low angular resolution observation. The LOSVDs are integrated through a circular aperture with a 
Gaussian weight function with a FWHM of 10, representing an angular resolution five times worse than the binary separation. The left panel is for q= \ and 
the right for 9 = 0. 1 . The lines of site are centred at.i: = ±10 andx = ±5, i.e. on both sides of the BBH. There is a clear shift in the line peaks in the q=\ case 
between the two sides, resulting from the rotation structure evident in the side-view maps (see Fig. 4). The shift is significantly smaller in the ^ = 0. 1 case. 



sphere, are on highly unstable orbits and disappear quickly from 
the system. Stars can enter or leave the Hill spheres only through 
an energy exchange with a fourth body, which may be very slow in 
low density cores. 



6 CONCLUSIONS 

Orbits in the restricted three body problem are notoriously complex 
(though some insight can be gained from stability maps). Here we 
exploit this property to derive the signature of a BBH on the nearby 
stellar kinematics, once unstable orbits are gone. The fraction of 
velocity phase space populated by stable orbits decreases inwards. 
The stars ejected by unstable orbits will leave behind a light defi- 
ciency, which was suggested to explain the core structure of mas- 
sive ellipticals. The remaining stars are on significantly anisotropic 
orbits, characterized by the following properties: 

(i) Tangential orbits dominate, mostly retrograde at r < 5. 
(ii) Increased CT, due to the elimination of low tangential veloc- 
ity orbits. 

(iii) A drop in a at r < 5, as most orbits become retrograde. 

These properties lead to a specific signature on the LOSVD mo- 
ments on scales as large as 5-lOiZhi which may be resolved in 
nearby galaxies. The detection of these kinematic features may in- 
dicate the presence of a BBH currently, or a relaxation time ago, 
beyond which the kinematic signature is erased. If the core struc- 
ture is a signature of a BBH phase in the past, some BBH kinematic 
signature may remain on the core scale as well. 
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